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Abstract. In Kelly et al. 0J, we presented new binary black-hole initial data adapted to 
puncture evolutions in numerical relativity. This data satisfies the constraint equations to 
2.5 post-Newtonian order, and contains a transverse-traceless "wavy" metric contribution, 
violating the standard assumption of conformal flatness. We report on progress in evolving 
this data with a modern moving-puncture implementation of the BSSN equations in several 
numerical codes. We discuss the effect of the new metric terms on junk radiation and 
continuity of physical radiation extracted. 
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1. Introduction 

The astounding success of numerical relativity in simulating the merger of comparable-mass 
black-hole binaries in recent years stemmed from a number of numerical approaches to 
initial data, evolution formulations, gauge conditions, and even grid structures. However, 
many active groups have converged on a simple combination of methods called the "moving 
puncture" prescription [3j . 

To initialise the numerical fields for a puncture evolution, most groups use the puncture 
prescription of Brandt & Briigmann 01 with Bowen-York extrinsic curvature [5]. In this 
scheme, the three-metric 7^ is conformally flat: 

Hj = ^Vij, (1) 
and the conformal factor ip must satisfy the Hamiltonian constraint 

A-0 + l -K ab K ab i,~ 7 = 0, (2) 

o 

where the conformal extrinsic curvature K^j already satisfies the momentum constraint for 
holes with arbitrary momentum and spin. The zero-momentum constraint can be solved 
exactly to yield the Brill-Lindquist solution for a pair of holes at a point of time-symmetry |[6l 

1{> = 1+ ZT7 + Jjp =7 ■ (3) 

2\x — X\\ 2\x — X2\ 

Here the are "bare" or "puncture" masses, residing at positions xa on the numerical 
grid. However, to solve Eq. in general requires numerical methods, and the infinities 
encountered at the puncture locations are problematical. Brandt & Briigmann's insight was 
that the divergent parts of ip could be formally factored out, leaving a well-behaved, simply 
connected sheet on which to solve their modified constraint. 

With the Brandt-Briigmann prescription, only a single elliptic equation has to be solved, 
and several fast solvers have been developed to perform this operation to extremely high 
precision. The physical mass Ma of each black hole after solution will be greater than its 
puncture mass tua- 

However, the restriction of this data is that it is, by construction, conformally flat. We 
know that the Kerr metric, the archetypal stationary solution of Einstein's vacuum equations, 
is not conformally flat unless it has vanishing spin. It would seem that requiring conformal 
flatness of an astrophysically realistic binary (which has significant orbital angular momentum 
by construction) is unrealistic. In practice, when evolving puncture binary data, we see 
early bursts of unphysical high-frequency radiation propagate through the domain; see Fig. [H 
reproduced from Q, as an example of this radiation at different initial separations. 

2. Post-Newtonian Metric in the ADM-TT Gauge 

In the 1970s, Ohta et al. jH HI [TOj derived conditions for the post-Newtonian metric of an 
iV-particle system in the trans verse-traceless ADM (ADM-TT) gauge. The structure of this 
metric was given to 2.5 PN order by Schafer |fl"Tfl: 

Jij = ippnVij + hJ^ > (4) 
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Figure 1. Strain waveforms for an equal-mass, zero-spin binary of total mass M, from 
three initial separations (/if indicates the resolution of the highest- refinement region in 
the numerical simulation supplying each waveform). Here the time axis is measured in units 
of the remnant mass ni{ < M of the post-merger hole, while the strain h is multiplied by the 
extraction radius R (again in units of mf). Note that the waveforms in the 9. 54M and 6.51M 
cases exhibit an initial "glitch". Adapted with permission from [7]. 



where the post-Newtonian conformal factor ip is expanded as 



if;- 
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(5) 



In this expression, 0( 2 ) alone yields the Brill-Lindquist potential term © for two stationary, 
nonspinning particles, while the leading corrections in <j)u\ depend explicitly on the separation 
of the particles, and on their momenta. Note that the three-metric 7^ is no longer conformally 
flat, due to the presence of a trans verse-traceless term hf^. This satisfies an outgoing wave 
condition: 

2- ~^ — 6{x-x A ) + -t 

IA=1 m A 4 
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TT 
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(6) 



The corresponding extrinsic curvature is derived from the conjugate post-Newtonian 
three-momentum: 
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Explicit expressions were given for the terms 0( 2 ), 0(4) and 7r^ by [1TT1I . who also suggested 
a "near zone" approximation for hf^, by splitting the retarded inverse d'Alembertian in © 
with an inverse Laplacian: 

" N 1 

^2 p u / ; i' - 1 



h 



TT 



^TT (JVZ) _|_ ^TT (remainder) 



TTkl 

ij 



A=l 



m A 



8{x 



(8) 
(9) 



+ 0(*;/c) 5 . 

The explicit form of this near-zone approximation was supplied by |[T2l|. 

Tichy et al. lfT3l adapted the ADM-TT-gauge 2.5PN results to puncture initial data 
for NR. They established that although the ADM-TT metric was not conformally flat, the 
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Figure 2. Descriptive plot of the portion of the trajectory (thin curve) of particle 1 relevant 
to the evaluation of h TT at field point P. hf„ T ^ NZ ^ and hf„ T (P resent *> ^ evaluated with the 
instantaneous positions and momenta at initial time t = 0. \h L J ( retarded ^ [ s evaluated using 
the instantaneous positions and momenta at the retarded time t — t r relative to P (found by 
travelling back along the light cone, indicated by a thick dashed line, until it intersects the 
particle's trajectory). Finally, hf^ ( m * e ™ a ') [ s evaluated by integrating along the trajectory 
segment from t = V to t = 0. 



behaviour of the metric near the black holes was dominated by the conformal factor, and so 
fixed-puncture evolution methods should work as for standard puncture data. 

A few years later Kelly et al. 0] completed the picture for nonspinning binaries by 
determining the "remainder" TT term, hj^ ( remamder *> i 2PN order. In the next section, we 
highlight some of the main properties of the complete solution. 



3. Summary of Global Properties of Solution 

While the work of [[TTi and |[T2l applies to general systems of particles, the "remainder" term 
presented in [IJ applies only to the simplified situation of a binary system (N = 2). In such a 
system, they determined that the structure of the remainder term divides into three segments, 
according to time of evaluation: 

yTT {remainder) _ yTT (present) _|_ yTT (retarded) _|_ yTT (interval) 

For each field point where hj^ is to be evaluated, the "present" term is evaluated using the 
particle positions and momenta at t — 0, the time at which the simulation will start. The 
"retarded" term is evaluated using positions and momenta at the retarded time of each source 
particle relative to the field point. Finally, the "interval" term is an integral over the particles' 
paths from the retarded time to the present. Figure [2]illustrates this division. 

Each of these segments, moreover, consists of a "kinetic" and a "potential" part, the 
former depending on the particles' momenta, the latter on their relative positions. 

In this expression, the present-time piece almost completely cancels the near-zone 
solution of [12J: the kinetic terms cancel exactly, while the slightly more involved potential 
terms are suppressed by three powers of the field distance R: 

u TT(NZ+present) G 2 m 1 m 2 r f I \ 

The retarded-time piece reduces to the well-known quadrupole solution for a nonspinning 
binary as r/R — > 0. The most involved term, the interval piece, is too difficult to do in 
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Figure 3. h TT (^ uH ) versus quadrupole for an equal-mass binary, evaluated along the orbital 
axis. Reproduced from JT]. 





Figure 4. Three-metric component j xx in the x-y plane with all terms for an equal-mass 
binary. 

generality, and must be integrated numerically. In Figure [3l we evaluate the full solution hf^ 
over time, along the orbital axis of an equal-mass system, assuming a simple inspiral. We see 
that the waveform is very close to the quadrupole solution. 

Figure |4] shows one of the three-metric components for the full solution, with the 
characteristic quadrupole swirl we expect from an inspiralling binary. 

4. Encoding the Binary's Past Inspiral 

m i , .< j j , ■ j . ■ • . i . •« , • ; TT (retarded) ( TT (interval) 

To evaluate the retarded-time and time-interval contributions, , , we 

need a model for the past history of the black holes. Initially, in HI, we employed a hybrid 



Post-Newtonian Initial Data with Waves: Progress in Evolution 



6 



procedure (here M = M x + M 2 is the total mass of the system, while 77 = M 1 M 2 /M 2 is the 
symmetric mass ratio): 

(i) obtain separation r as a function of orbital frequency to 2PN order by (numerically) 
inverting the "puncture adapted" relation from [fl~4|| : 

MVL 



64M 3 r 3 1 M A r] 1 M 5 (-5rj + 8ri 2 ) _ 
\\ (M + 2r) 6 + c 1 ^ r+ c 1 8^ ' ( } 

(ii) obtain transverse momentum p as a function of f2 to 2PN order from Schafer & Wex 
E9: 

6 72 

(iii) obtain orbital phase $ (and hence frequency fi = d$/ c?t) as a function of time to 2PN 
order using the explicit relation from lfT6l : 

05 /8 /3715 55 \ /8 _ Stt /4 1 
V8064 96 7 4 J ' 

where = 77 (t c — t)/5 M, and i c is the (nominal) merger time. 



77 



This method has several drawbacks. For one thing, it is quite limited in post-Newtonian 
order. For another, the components are in inconsistent gauges. Finally, we have no 
prescription for an instantaneous radial momentum, necessary for low-eccentricity inspiral. 

A conceptually simpler approach is to get all the needed information from a single source. 
Following recent practice in initial parameters for numerical evolutions of punctures [17J, we 
can evolve the binary system inspiral through Hamiltonian evolution of the PN equations of 
motion. This has been shown by Husa et al. [17] to result in extremely low eccentricity, 
at least for nonspinning data. For more generic spinning binaries, the situation is more 
complicated, but promising; see, for instance, |[T8l . 

Although we need much more information for the new data, the Hamiltonian evolution 
method should perform just as well as for simple punctures: the puncture positions and 
momenta required are all from earlier times, and hence larger separations with lower 
velocities, where post-Newtonian methods are guaranteed to work. The only drawback is 
that potentially a lot of data must be stored about the past history of the binary to calculate 
the retarded and interval terms: and the larger the numerical grid, the further back in time we 
must reach. 



5. Numerics 



In the geometric units (G — c — 1) commonly used for vacuum numerical relativity, time, 
length, and mass can be scaled by a single number. For this purpose we use M = Mi + M 2 , 
the total mass of the binary system. 

The numerical implementation of the wavy PN data has taken place in three independent 
codes, the Cactus-based LazEv code 0HI21I2Q1, the BAM code [EH [22l [23l , and the Paramesh- 
based Hahndol code [|24ll25l . Some of the code was auto-generated using aMathemat ica 
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Figure 5. Left panel: Hamiltonian constraint in the orbital (x-y) plane, with (dashed) and 
without (solid) h^J i mterval ) Right panel: Hamiltonian constraint in the same plane with full 
h TT (dashed) and with no h TT (solid). 

script supplied by Gerhard Schafer. For simplicity, the time-derivatives of hf^ appearing 
in the 2.5PN extrinsic curvature calculation CO) are carried out using simple second-order- 
accurate centred differencing (this is easily accurate enough, as the time spacing used is much 
smaller than the spatial discretisation of the numerical grid). 

To calculate the two retarded times for each field point (one for each black hole), we use 
a Newton solver. Interval terms are then calculated by integrating from these retarded times 
to the present using Romberg integration. 

Before discussing the evolution of the data, we note two numerical properties of the 
initial data. The first is that, as expected, the Hamiltonian constraint violation for the complete 
solution is better than for a partial solution: the left panel of Fig. [5] demonstrates this with the 
Hamiltonian constraint evaluated for the complete solution with and without the "interval" 
term. 

On the other hand, it seems that leaving out the h TT terms entirely yields even lower 
constraint violations, as seen in the right panel of Fig. [5l The reason for this is unclear, but 
may have to do with the greatly increased number of numerical evaluations in the h TT terms. 
Given the observed behaviour of the black-hole (horizon) masses described in the next section, 
it seems that the bulk of the constraint violation "falls back" into the hole, with little escaping 
to the field zone to pollute the waveforms. We discuss the residual constraint violation further 
in Sec. El 

6. Early Evolution Results 

Initial evolutions have been carried out for equal-mass nonspinning binaries, with initial 
separations between 6M and 10M. All evolutions show a combination of desirable and 
undesirable effects, which we illustrate with figures from lOM-separation evolutions at low 
central resolution (3M/128 near the punctures) with the Goddard Hahndol code. We will 
often compare with the results of a more traditional moving puncture evolution of Bowen- 
York puncture data, with initial separation of 1 1 M, using parameters from ifTTl . For this data, 
we solved the Hamiltonian constraint numerically using Marcus Ansorg's TwoPunctures 
spectral code E6l . In both cases, we tracked apparent horizons before and after merger with 
Jonathan Thornburg's AHFinderDirect code ll27ll . 
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Figure 6. Coordinate separation of traditional "solved" data (solid) and wavy data (dashed). 
The wavy data exhibits almost 10% eccentricity initially. 



6.1. Eccentricity and Horizon Masses 

The first thing we note is the presence of strong eccentricity in the puncture tracks of the holes 
- see Fig. [6] This eccentricity appears to be around 10%, far greater than that of the traditional 
evolution, and persists until around 100M before merger. 

A related phenomenon may be seen in the evolution of apparent horizon masses for the 
pre-merger binary - see Fig. [7] The apparent horizon mass Mah of a hole is related to the area 
of the apparent horizon. We locate the former numerically using the AHFinderDirect 
code ll27l . We calculate the proper area Aak of this 2-surface, and derive from this the 



irreducible mass M irr = J Aah/IQtt. The horizon mass is then related to M irr by inverting 
the Christodoulou formula 1281: 



2Ml 



1-y/T 



(15) 



where a = \S\/M 2 is the dimensionless spin of the hole. Note that for zero-spin holes, the 
horizon mass Mah is identical to the irreducible mass M irr . 

In Fig. |7] we see that while the standard horizon masses quickly settle down to a stable 
value, varying insignificantly over the 1100M of pre-merger evolution, the unsolved wavy 
data masses decay at a much slower rate, and with a periodic saw-tooth feature over time. 
If this is not merely an artifact of the horizon- finding algorithm or implementation, then the 
holes are losing mass steadily over the course of the inspiral. As discussed in [T30l . this could 
be associated with the nature of the residual constraint violation on the initial time-slice. This 
would lead to a considerable ambiguity in what the "correct" horizon mass is. The initial 
momenta (as well as the entire past history in the h TT terms) depend sensitively on the hole's 
mass; an incorrect mass might result in considerable eccentricity. 
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Figure 7. Pre-merger AH mass for black hole 1 for traditional data (solid) and wavy data 
(dashed). Over the course of the inspiral, the wavy data horizon mass drops by around 1%. 
Variation in the mass for the traditional data is negligible by comparison. 



6.2. Waveforms 

The main quantity of interest from these evolutions must be the gravitational radiation 
extracted. In the left panel of Fig. [8] we show the real part of the dominant (£ — 2,m = 2) 
waveform mode Rip4, as extracted on the coordinate sphere R = 45M. We note that indeed 
the physical waveform is present from t = 0, with less junk radiation than is present in the 
standard solved data. 

Unfortunately a poor choice in the structure of the numerical grid in the radiation zone 
led to noisy extraction for both traditional and new data. Note that a high-accuracy waveform 
would be extracted at (or extrapolated to) R — > oo; however, this was not possible with the 
numerical grid used for these initial evolutions. Coupled with the relatively low resolution 
used, we estimate waveform amplitude errors of up to ~ 8% at the peak; thus our numerical 
results are qualitative in nature at this point. Higher-accuracy extraction and analysis methods 
will be appropriate in future evolutions, when the mass and eccentricity issues described above 
have been resolved. 

6.3. Final State 

The final state of the post-merger black hole, which we analyse using the AHFinderDirect 
code ll27ll . is qualitatively similar to that of the standard solved data. The final spin, as 
estimated by values of the Coulomb scalar on the surface of the apparent horizon ||29l , is 
around a = S z /M 2 w 0.68625. Using (TT5T) . we estimate the total horizon mass Mah of the 
remnant hole. Figure [9] shows a and A/ah f° r the standard and wavy merger remnants. 



Post-Newtonian Initial Data with Waves: Progress in Evolution 



10 




100 200 300 200 400 600 



Figure 8. (2,2) mode of Rip4 from t = for standard puncture evolution (solid) and wavy 
data (dashed), as extracted at R = 45M. The left panel shows the real part of the mode, 
which already has a non-zero value at t = 0, and a smaller non-physical pulse around 
t = 45M, compared with the traditional moving-puncture waveform. The right panel shows 
the amplitude of complex mode. The high eccentricity visible in the (gauge-dependent) 
puncture tracks has a clear effect here also, compared with the low-eccentricity traditional 
evolution. 
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Figure 9. Left panel: post-merger z-spin a = S z /M 2 for standard data (solid) and wavy 
data (dashed). Right panel: Post-merger horizon mass for standard data (solid) and wavy data 
(dashed), using Eq. (JT3J, and assuming a = 0.68625. 



7. Open Questions & Future Work 

From these early evolutions, we see that the the new wavy PN data appears to achieve at least 
some of our goals: it does evolve stably in the moving puncture recipe, without any special 
tweaks in gauge or evolution equations. Moreover, the early waveforms do indeed show more 
reasonable start-up behaviour than those of a standard puncture evolution, with a physically 
reasonable non-zero value and a diminished nonphysical pulse. 

Nevertheless, some aspects of the numerics are not satisfactory, most notably the high 
orbital eccentricity and the slowly settling horizon masses of the pre-merger holes. These two 
phenomena may be simply related, as our experience in working with standard puncture data 
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has shown us that small errors in tuning masses will lead to eccentricity. It is conceivable that 
the residual Hamiltonian constraint violation of the wavy data causes this mass defect; similar 
effects have been studied recently in ll30ll . 

Resolving the mass issue may necessitate the introduction of a numerical elliptic solver 
to remove the residual constraint violation. To do this completely is not straightforward, as 
the ADM-TT gauge used here will change the form of the Hamiltonian constraint; moreover, 
the momentum constraint will in general need to be solved too. 

Looking beyond these issues, we may consider the introduction of spin to our data. 
Though explicit post-Newtonian solutions of the constraint equations with spin are not yet 
available, we note that the leading-order momentum contributions to the conformal extrinsic 
curvature are just the Bowen-York momentum terms. It is conceivable that the leading-order 
spin contributions are given by the corresponding Bowen-York terms also. 

Finally, we note that we have not addressed the initial conditions of the lapse function 
and shift vector. It is known that in standard puncture evolutions with the popular "1+log" 
slicing conditions, these settle down at late times to a "trumpet" form OTl l32ll . Until these 
late-time forms can be incorporated into the full initial data, we cannot expect to eliminate 
gauge pulses in our waveforms. 
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